Role of the rhizosphere bacterial community in assisting phytoremediation in a lead-zinc area

Heavy metals (HMs) contamination and vegetation destruction in the mining area caused by mining activities are severely increasing. It is urgent to restore vegetation and stabilize HMs. In this study, we compared the ability of HMs phytoextraction/phytostabilization of three dominant plants, including Artemisia argyi (LA), Miscanthus floridulus (LM), and Boehmeria nivea (LZ) in a lead-zinc mining area in Huayuan County (China). We also explored the role of the rhizosphere bacterial community in assisting phytoremediation using 16S rRNA sequencing technology. Bioconcentration factor (BCF) and translocation factor (TF) analysis showed that LA preferred accumulating Cd, LZ preferred accumulating Cr and Sb, and LM preferred accumulating Cr and Ni. Significant (p < 0.05) differences were found among the rhizosphere soil microbial communities of these three plants. The key genera of LA were Truepera and Anderseniella, that of LM were Paracoccus and Erythrobacter, and of LZ was Novosphingobium. Correlation analysis showed some rhizosphere bacterial taxa (e.g., Actinomarinicola, Bacillariophyta and Oscillochloris) affected some soil physicochemical parameters (e.g., organic matter and pH) of the rhizosphere soil and enhanced the TF of metals. Functional prediction analysis of soil bacterial community showed that the relative abundances of genes related to the synthesis of some proteins (e.g., manganese/zinc-transporting P-type ATPase C, nickel transport protein and 1-aminocyclopropane-1-carboxylate deaminase) was positively correlated with the phytoextraction/phytostabilization capacity of plants for heavy metals. This study provided theoretical guidance on selecting appropriate plants for different metal remediation applications. We also found some rhizosphere bacteria might enhance the phytoremediation of multi-metals, which could provide a reference for subsequent research.


Introduction
Due to the continuous mining activities, heavy metals (HMs) contamination and vegetation destruction are severely increasing. It is a severe threat to ecological biodiversity and human health. Therefore, soil restoration and governance need urgently be solved (Xiang et al., 2021;Zerizghi et al., 2022). The primary remediation methods for metal-contaminated soils include physical, chemical and biological strategies.
In general, physical and chemical remediation methods are costly, unfriendly to the environment, and can easily lead to secondary pollution, so that many researchers now focus on bioremediation strategies (Gonzaĺez Henao and Ghneim-Herrera, 2021;Azhar et al., 2022). Phytoremediation has been widely used due to its advantages of being cost-effective, economical, eco-friendly, and sustainable (Ahemad, 2019;Shrivastava et al., 2019). Currently, As hyperaccumulators such as Pteris vittata and Pteris cretica, Zn hyperaccumulator Sedum alfredii, Cd hyperaccumulator Viola baoshanensis, and Mn hyperaccumulator Phytolacca acinose have been found in China for the phytoremediation of mine tailings (Niu et al., 2021). The effectiveness of phytoremediation depends on the chemical and physical properties of the plant, the bioavailability of metals in the soil, and the ability of rhizosphere soil microorganisms to absorb, transfer and detoxify metals (Subpiramaniyam, 2021). The removal of HMs by microorganisms has the advantages of easy use, low cost, large adsorption capacity, and high efficiency. Among these, bacteria, fungi and algae are widely used (Yin et al., 2019). In general, heavy metal ions can be adsorbed and combined by some functional groups of bacterial polysaccharide mucus layer, e.g., carboxyl, amino, phosphate, and sulfate (Yue et al., 2015). The adsorbed heavy metal ions can enter the microbial cells by metal-related enzymes or proteins, changing the redox state of heavy metal ions and thereby reducing their toxicity.
In practical applications, plants and microorganisms for bioremediation are more efficient (Vaid et al., 2022). Usually, microorganisms can improve plant extraction by increasing the availability of HMs in plants and increasing plant biomass (Wood et al., 2016;Mishra et al., 2017). Previous studies have shown that Patescibacteria increased the availability of HMs in the rhizosphere and promoted the remediation of heavy metalcontaminated soil by Sedum alfredii . Moreover, the plant preference for some HMs is influenced by the rhizosphere ecological characteristics, including the specific hormones, root exudates, soil nutrients, soil properties and rhizosphere soil microbes (Bais et al., 2006). The phytoremediation efficiency was regulated by soil enzyme activity and beneficial rhizosphere-associated microorganisms Trifolium repens L. (Lin et al., 2021). Therefore, the combined plant-microbes method can improve the heavy metal resistance of plants and achieve an ideal remediation effect. In recent years, it has also been found the limited efficiency of phytoremediation with a single plant (Buscaroli et al., 2017;Rizwan et al., 2019;Hosseinniaee et al., 2022) and the co-planting pattern of complementary plants for metals enrichment may be more efficient. The co-plantation of Solanum nigrum with Quercus nuttallii or Quercus pagoda effectively improved the enrichment of cadmium (>40%) and zinc (>30%) (Qu et al., 2021).
The environmental problem is severe in the lead-zinc mining area in Huayuan County, Hunan Province, China. This area is a super large ore deposit with many tailings ponds. Wastewater and ore sand was discharged through tunnels and village ditches, accumulating in the surrounding crops and villagers. Barbaric mining and poor management caused environmental pollution, damaged mine landforms, reduced vegetation coverage, soil erosion, and other problems. Therefore, we aim to explore the relationships among metals uptake of the dominant plants, rhizosphere bacterial community and the soil environment under the long-term HMs stress, analyze the ecological characteristics of the rhizosphere (physicochemical properties, soil enzyme activities and bacterial community structure), and the major microbial groups and metabolisms responsible for the oxidoreduction/ detoxification/tolerance of HMs in the rhizosphere and absorption/translocation of HMs of plants.

Sample collection and pretreatment
The Pb-Zn mining area was located in Huayuan County, Hunan Province (China), with a subtropical monsoon climate with an average temperature of 16.0°C and an average rainfall of 1363.8 mm. Environmental conditions and vegetation patterns were studied through field sampling. The vegetation pattern and distribution at different study sites were analyzed by quadratic analysis to identify dominant plant species, including Artemisia argyi (LA, 109°22′31.31″ E and 28°31′59.47″ N), Boehmeria nivea (LZ, 109°22′31.31″ E and 28°31′41.25″ N) and Miscanthus floridulus (LM, 109°21′41.12″ E and 28°30′46.47″ N) ( Figure S1). There were eight replicates for each plant and the corresponding rhizospheric soil. Soil samples were collected from rhizosphere soil (5-20 cm). All soil and plant samples were sealed into clean polyethylene bags and transported to the laboratory at a low temperature (<5°C ). Fresh plants were rinsed with tap water and deionized water. Then the plant samples were dried at a high temperature of 105°C for 30 minutes and then dried at 65°C to constant weight. The dried plant tissue was then separated into roots and shoots, ground to a fine powder, and placed in polyethylene bags for further analysis. Soil samples were airdried, ground to a particle size of less than 0.147 mm , and stored in polyethylene bags until analysis. The remained soil samples were stored at -80°C for microbial analysis.
Soil pH was measured by a pH meter (PHS-25, Shanghai Yidi Scientific Instrument Co., Ltd., China). 5.0 g soil was dried at 65°C to constant weight to measure the moisture content (MC). Soil-available phosphorus (AP), soil-available potassium (AK), soil acid phosphatase (S.ACP), soil catalase (S.CAT) and soil peroxidase (S.POD) were measured with the test kit (Shanghai zcibio technology co., Ltd. China). Analysis of soil organic matter (OM) was measured as described by Walz et al. (Walz et al., 2017). Total nitrogen was determined using the Kjeldahl method (Bremner, 2009). Soil HMs were determined regarding the method used by Sungur et al. (Sungur et al., 2020).

DNA extraction and Illumina MiSeq sequencing
DNA extraction and Illumina MiSeq sequencing of soil microorganisms were performed by Shenzhen E-Gene Technology Co., Ltd., using a 16S rRNA high throughput sequencing technique. Briefly, PCR was performed using primer pair 515 F (50-GTGCCAGCMGCCGCGGTAA-30) and 806 R (50-GGACTACHVGGGTWTCTAAT-30). PCR reaction conditions were as follows: 3 min of denaturation at 95°C , 27 cycles of 30 s at 95°C , 30 s for annealing at 55°C , and 45 s for elongation at 72°C , and a final extension at 72°C for 10 min. The raw data were deposited in the Sequence Read Archive of the National Center for Biotechnology Information database (accession number: PRJNA904575).

Molecular ecological network 2.4.1 Network construction
In this study, the relative abundances of bacterial OTUs in the soil were used to construct phylogenetic molecular ecological networks (PMEN) using a network approach based on random matrix theory (RMT). This method can automatically identify keystone OTUs and determine the topological properties of the network. The detailed operation methods were referred to the previous study by Deng et al. ( 2012). In our study, the OTUs detected in equal or more than 3 of the 8 biological replicates were kept for network construction. The default similarity threshold is 0.98. Finally, the network interactions were visualized using Cytoscape 3.9.1 and Gephi 0.9.5.

Module detection and keystone microorganism identification
Molecular ecological networks consist of many OTUs/genes (nodes) and interactions (links) (Dunne et al., 2008). A module is a group of nodes that have similar functions and effects. Maximum modularity can separate the OTU into multiple dense subgraphs (Tao et al., 2018). Within a module, the role of a node is characterized by its intra-module connectivity (Zi) and intermodule connectivity (Pi). Based on Zi and Pi, nodes in the network could be classified into four different roles: peripheral nodes (Zi ≤ 2.5 and Pi ≤ 0.62, nodes have only a few links with other nodes within their modules); connectors (Zi ≤ 2.5 and Pi > 0.62, nodes highly connected to several modules); module hubs (Zi > 2.5 and Pi ≤ 0.62, nodes highly connected to many nodes within their own modules); and network hubs (Zi > 2.5 and Pi > 0.62, nodes act as both module hubs and connectors (Olesen et al., 2007). Nodes playing the roles of connectors, module hubs, and network hubs might act as keystone taxa to maintain the network structure and function (Faust and Raes, 2012).

Functional profiling
Before functional gene prediction using PICRUSt (phylogenetic investigation of communities by reconstruction of unobserved states) described by Langille et al. (2013), the detected OTUs were reclassified using the GREENGENES reference database. Subsequently, PICRUSt used 16S rRNA genes to infer metagenome gene functional content from phylogenetic information. The predictions were precalculated for genes in databases, including the Kyoto Encyclopaedia of Genes and Genomes (KEGG). The input data were first normalized by copy number by dividing each OTU by the known 16S copy number abundance before metagenome predictions and subsequent collapse into different functional pathways. The output of PICRUSt consisted of a table of functional gene counts as KEGG orthologs (KOs). The Nearest Sequenced Taxon Index (NSTI) value was used to validate the reliability of predicted metagenomes and functional pathways.

Statistical analysis
A completely randomized experimental design with eight replicates per treatment was used. Statistical analysis was conducted using analysis of variance (Markowitz et al., 2014). Multiple comparison analyses among treatments were performed using IBM SPSS Statistics 26 with Tukey's test. The translocation factors (TF) and bioconcentration factors (BCF) of three plants were calculated to evaluate the restoration potential of plants. TF evaluates the ability of plants to transport and enrich HMs from underground to aboveground. BCF reflects the ability of plants to absorb and transfer HMs into plant tissues.
BCF and TF were calculated as follows: BCF= C plant/C soil; TF= C shoot/C root.
The Shannon and Simpson diversity index analysis showed species abundance and evenness. The Chord diagram of the species relationship map was used to analyze the species composition of each sample in the microbial community. A Venn diagram was used to show the general relationship between the different treatments. It was commonly used to visualize common or unique information from multiple samples. Ternary phase diagrams showed the relative abundance of species in different groups. Principal Component Analysis (PCA) and Non-metric multidimensional scaling (NMDS) analysis were used to reflect the distances of the bacterial communities among samples. LEfse analysis (LDA Effect Size analysis) was conducted to find the species with significant differences in abundance between groups. Correlation analysis was conducted to show the Spearman correlation coefficient among environmental variables and the Mantel correlation coefficient among microbial species, predicting functional genes, environmental variables and BCF/TF. All the above mapping analysis was conducted using Origin 64 and Rstudio 4.2.0.

The contents of HMs in the shoots and roots.
The contents of HMs in the shoots and roots of plants were shown in Table 1. For the shoot, the content of Pb in the LA group reached 144.43 mg/kg, significantly (p < 0.05) higher than that in the LZ and LM groups. The content of Zn (635.81 mg/ kg), Cu (81.29 mg/kg), Fe (5529.07 mg/kg) and Mn (501.26 mg/ kg) in the LZ group was significantly (p < 0.05) higher than that in the LA group. The content of Cr (654.82 mg/kg) and Ni (90.8676 mg/kg) in the LM group was significantly (p < 0.05) higher than that in the LA and LZ groups. Except for Mn, the contents of other HMs in the roots of LA were significantly (p < 0.05) higher than those of LZ and LM in the root. Among them,

Soil physicochemical/biochemical properties
The physicochemical properties in the rhizosphere soil were shown in Table 2. The MC (2.75%-15.88%), pH (4.89-6.37), OM (3.47-21.55) and TN (0.36-2.29) were different in rhizosphere soils of different plants, and they were significantly (p < 0.05) higher in the LZ group than those of LM. The soil AK contents showed no significant (p > 0.05) difference among the three groups, but the soil AP content in the LZ and LM groups was significantly (p < 0.05) higher than that of LA by 18.77% and 22.53%, respectively. The activity of S-ACP in the LA and LZ groups were significantly (p < 0.05) higher than that in the LM group, but the activities of S-CAT and S-POD showed no significant (p > 0.05) differences.
In addition, the contents of As (11.17 mg/kg), Cd (2.53 mg/ kg), Cr (3.83 mg/kg), Cu (29.32 mg/kg), Fe (14148.06 mg/kg), Mn (761.88 mg/kg), and Ni (21.69 mg/kg) in the LA group were significantly (p < 0.05) lower than those in the LZ or LM group. The Pb content in the LZ group was 392.18 mg/kg, which was significantly (p < 0.05) lower than that of LA and LM.

Translocation factor and Bioconcentration factor
The results of BCF and TF were shown in Table 3. The BCFs of Cd (2.31), Cr (105.40), Cu (2.93) and Ni (1.66) were greater than 1 in the LA group, while only the BCFs of Cr (11.04) and Sb (1.74) in the LZ group and the BCF of Cr (16.36) and Sb (3.08) in the LM group were greater than 1. In addition, LZ also has strong TF (TF > 1) for all other metals, LM has strong TF (TF > 1) for all metals except Mn and Sb, and LA has strong TF (TF > 1) for Cd and Zn. The TFs of Cr, Cu, and Ni in LA also reached 0.84, 0.71, and 0.83, which were relatively higher than the TFs of other HMs in the same plant species. It was found that the TF of LZ for Cd, Cu, Mn, Ni, Pb, Sb, and Zn were higher than those of LA and LM, while the TF of LM for As, and Cr, was the strongest.

Overview of bacterial community in the rhizosphere soil
Simpson Index showed no significant (p > 0.05) differences among the three groups, but the Simpson Index in the LM group reached 65.51, which was higher than the other two groups. The Shannon Index showed that the microbial diversity of LA, LM, and LZ were 4.11, 4.18, and 4.44, respectively. Among them, LM had the most abundant diversity and was significantly higher than LA ( Figure 1A). According to the Venn diagram, there were 26 duplicated OTUs in 864 OTUs ( Figure 1C). The abundance chord chart showed that soil bacterial communities mainly consisted of Armatimonadetes_gp4, Gp3, and Thermanaerothrix, and the dominant genera were different among three groups ( Figure 1B): they were Ornatilinea (4.47%), Gp16 (3.63%), and Saccharibacteria_genera_incertae_sedis (3.60%) in the LA group; they were Armatimonadetes_gp4 (6.92%), Gp3 (4.93%), and Thermanaerothrix (3.43%) in the LM group; They were Gp3 (4.66%), Bacillariophyta (4.30%), and Iphinoe (3.50%) in the LZ group. A ternary plot analysis performed on the phylum level showed that Actinobacteria, Cyanobacteria/Chloroplast, Chloroflexi, Firmicutes, and Proteobacteria were the dominant phyla ( Figure 1D). Armatimonadetes dominated LM, and Planctomycetes dominated LZ. Results of PCA showed that points of samples in the three groups separated from each other group ( Figure 1E), and similar results were obtained by dissimilarity tests, suggesting that significant (p < 0.05) differences occurred in the rhizosphere bacterial community structure of different plant (Table S1).
Further LEfse analysis at the genus level was shown in Figure 2. The cladogram (Figure 2A) compares the differences of genera between three plants, and LDA analysis ( Figure 2B (E) Principal Component Analysis. In the ternary phase diagram, different points represent different species, the size of the points represents the average abundance of the species in different groups, the gray points represent the non-enrichment, and the points of each color respectively represent the enriched groups. shows the main genera with significant (p < 0.05) differences. The relative abundances of genera (e.g., Ornatilinea, Gp16, Enteractinococcus) were highest in the LA group; The relative abundances of genera (e.g., Armatimonadetes_gp4, Aggregatilinea, Paracoccus) were highest in the LM group; The relative abundances of genera (e.g., Bacillariophyta, Iphinoe, Gemmatirosa) were highest in the LZ group.

Co-occurrence networks analysis
Microbial population data were analyzed using the RMT-based network approach to discern the ecological network. Three networks were constructed based on the 16S rRNA gene sequencing data of three plans, respectively ( Figure 3A). Major topological properties of the empirical MENs of microbial communities in the eight groups are shown in Table S2. With the same threshold (0.980), their correlation values were more than 0.580, indicating that the degree of distributions in the constructed molecular ecological networks fits the power law model well. There were 134 links and 67 nodes in the LA group, of which the positive links accounted for 70.07%, and the negative links accounted for 29.93%. There were 137 links and 80 nodes in the LM group, of which the positive links accounted for 53.71%, and the negative links accounted for 46.29%. There were 175 links and 88 nodes in the LZ group, of which the positive links accounted for 63.52% and the negative links accounted for 36.48% ( Figure 3A). The network indexes showed that the microbial network of LZ was the most complex, but the number of modules in the LM group was higher than that of LA and LZ (Table S2). Results of ZP analysis showed that there were 6, 7 and 27 key OTUs in the LA, LM and LZ groups ( Figure 3B). Sub-network analysis of key OTUs revealed that OTU167 (Geothermobacterium, 4 connections, 4 positive correlations) and OTU89 (Clostridium XlVa, 5 connections, 4 positive correlations) were the core OTUs in the LA group. OTU77 (Paracoccus, 8 connections, 7 positive correlations) and OTU12 (Algisphaera, 4 connections, 4 positive correlations) were the core OTUs in the LM group. The core OTUs in the LZ group were OTU57 (Iphinoe, 19 connections, 16 positive correlations), OTU61 (Paracoccus, 18 connections, 14 positive correlations), OTU268 (Pseudoscardovia, 6 connections, 6 positive correlations), OTU111 (Limosilactobacillus, 7 connections, 6 positive correlations) ( Figure 3C). Venn diagram analysis of essential microorganisms showed that 10 key OTUs were shared by LA, LM, and LZ ( Figure 3D), which belonged to the genera, including Armatimonadetes_gp4, Gp3, Algisphaera, Clostridium XlVa, Robinsoniella, Taonella, Verrucosispora, Erythrobacter, and Thermanaerothrix.

Correlation analysis
According to the correlation analysis ( Figure 4A), it was shown that the key microorganisms were mainly related to the BCF of Cu, Cr, Cd, Fe, Ni and Zn, and the TF of As, Cd, Cr, Cu, Mn, Ni, Pb and Zn. S.ACP, S.CAT, pH, MC, OM, AP and TN were also significantly (p < 0.05) correlated to the key microorganisms. The diversity of the bacterial community was mainly related to the TF of Ni and Sb, as well as MC and AP. Bacterial community structures were significantly (p < 0.05) correlated to the TF of Ni, the activity of S.ACP, soil pH, and the contents of MC. Most HMs were mainly related to the contents of AP and MC in soil ( Figure 4A). TF_Cd was significantly (p < 0.05) positively correlated to the contents of OM/TN, soil pH, and the activity of S.ACP. TF_Fe was significantly (p < 0.05) positively correlated with TN, BCF_Mn, and BCF_Ni were significantly (p < 0.05) negatively correlated with the content of OM (p < 0.05). There was no significant (p > 0.05) correlation between BCF_Sb and soil physicochemical properties. Similarly, there was no significant (p > 0.05) difference between the TF and BCF of HMs and the soil contents of AK. Further correlation tests explored the relationships among microorganisms, HMs,

FIGURE 2
Linear discriminant analysis effect size (LEfSe). (A) The evolutionary branching diagram, and (B) the distribution histogram. Xiao et al. 10.3389/fpls.2022.1106985 Frontiers in Plant Science frontiersin.org and physiochemical properties in soil. It was found that the essential microorganisms were mainly negatively correlated with S.CAT, positively correlated with the TF of HMs, and had no significant (p > 0.05) correlation with S.POD. Interestingly, the Bacteroidetes was significantly (p < 0.05) negatively associated with soil nutrients. Furthermore, most microbiotas were negatively correlated with the BCF of HMs. The correlations between microorganisms and environmental factors/BCF/TF were shown in Figure 4B. Some microorganisms (e.g., Actinomarinicola, Bacillariophyta and Oscillochloris) were positively related to organic matter, pH and the TF of metals, while were negatively related to the BCF of metals.

Analysis of predictive functional genes
NMDS analysis showed significant inter-group differences in the predictive functional genes contained in soil microorganisms from different plant rhizosphere, while intra-group differences were slight ( Figure 5A). There were significant (p < 0.05) differences in the relative abundances of some functional genes in different microbial communities (Table S3). The relative abundances of the genes about 1-aminocyclopropane-1carboxylate deaminase (ACC deaminase), tryptophan synthase alpha/beta chain, nickel transport protein and manganese/Zinc- transporting P-type ATPase C were the most abundant in LA. The relative abundances of the genes about nonribosomal peptide synthetases were the highest in LZ (Table 4). Furthermore, some genes were further selected to make fitting curves with the contents of HMs in plant root tissues ( Figure 5B). The relative abundance of the gene about manganese/zinc-transporting P-type ATPase C was significantly (p < 0.01) positively correlated with the contents of Cd, Zn, and Pb, and the relative abundance of the gene about nickel transport protein was significantly (p < 0.01) positively correlated with the contents of Ni.

Discussion
In recent decades, our demand for metal resources has been increasing, which accelerated the mining of metal mines. The continuous expansion of mining areas led to more and more arable land facing environmental pollution (Diaz-Morales et al., 2021;Li et al., 2022). The increase in HMs will stimulate the production of a large number of reactive oxygen species, which will seriously impact the quality of crops and vegetables (Clemens and Ma, 2016) and lead to human health risks (Xiang et al., 2021). Phytoremediation is a sustainable approach to remediating contaminated sites (Hasnaoui et al., 2020). Thus, the study of the interaction between the hyperaccumulators and microorganisms can help formulate effective remediation strategies for mining areas (Xuan et al., 2021). In this study, we investigated the rhizosphere ecological characteristics of the dominant plants around the mining area and further explored their potential for multi-metal(loid)s phytoremediation. It was found that As, Cd, Cu, Pb and Zn have reached potentially hazardous levels (Table 1). We found that Artemisia argyi, Miscanthus floridulus, Boehmeria nivea are the dominant enrichment plants in this area, and the application of dominant plants combined with microbial communities for the bioremediation of soil pollution was explored ( Figure 6).
Microorganisms can directly or indirectly change the availability of metals and promote plant growth and absorption of HMs by regulating soil physicochemical properties or secretion of secondary metabolites. Our study found that Cyanobacteria/Chloroplast, Chloroflexi and Acidobacteria were key phyla for the plant to accumulate HMs. Different plants had different dominant rhizosphere flora, which affected the remediation of different HMs. The dominant phylum Cyanobacteria/Chloroplast in LZ was significantly (p<0.05) positively correlated with TF of Cd, Cu, Mn, Pb and Zn ( Figure S2). However, few studies reported the responses of Cyanobacteria/Chloroplast and their roles in phytoremediation under HMs stress, and it might need further study to verify. Chloroflexi was the dominant phylum in the LA group ( Figure 1B), significantly (p < 0.05) higher than in LM and LZ ( Figure 2B). Studies have found that Chloroflexi had an obvious advantage in polluted soil (Hemmat-Jou et al., 2018;Koner et al., 2022), which could rapidly adapt to heavy metal stress and affect the content of organic matter to control the practical availability of Cr and Pb for plants (Tang et al., 2019). It is consistent with this study, where Chloroflexi was also significantly (p < 0.05) positively correlated with the organic matter ( Figure S2). Acidobacteria was the dominant phylum in the LM group ( Figure 1B) and was significantly (p < 0.05) positively correlated with BCF of Cr ( Figure S2). Previous research found that the relative abundance of Acidobacteria was inversely related to pH (Debnath et al., 2016). Sun et al. (2022) found that the relative abundance of Acidobacteria was positively correlated to the contents of total Cr and available Cr, and Acidobacteria could cut down soil pH. The key genera of LA

A B
FIGURE 5 NMDS analysis of the structure in predicting genes (A) and the correlations between some genes and the contents of heavy metals in plant roots (B).
were Truepera and Anderseniella, that of LM were Paracoccus and Erythrobacter, and of LZ was Novosphingobium. The oxidation-reduction and nitrogen fixation activities of Truepera allow it to support plant growth (Zhou et al., 2022). Anderseniella can eliminate metabolic wastes, heavy metals, and aromatic chemicals (Karimi et al., 2019). Paracoccus can produce acid, alter the chemical and physical parameters in the soil and encourage plants to absorb heavy metals (Carvalho et al., 2018).
Erythrobacter is an iron metabolism bacterium that can secrete iron carriers and promote plant growth (Li L. et al., 2020). Novosphingobium has strong antioxidant activity and can slow down the toxic effect of heavy metals on plants (Petruk et al., 2018). Except for the dominant flora of different plants, the study found that they shared the same vital microorganisms ( Figure 3D). These microorganisms (e.g., Proteobacteria, Acidobacteria and Firmicutes) could assist plants in accumulating HMs. Studies have found that in the harsh tailings environment, the colonies mentioned above belonged to the dominant flora and could benefit the expression of metal resistance genes (Guo et al., 2017;Jiang et al., 2021;Koner et al., 2022). HMs in the soil can affect the growth of microorganisms through protein denaturation, cell membrane damage and inhibition of RNA expression and metabolism Duan et al., 2022). It has also been studied that HMs stress enhanced the functional fitness of endophytic bacterial communities (Yao et al., 2022). Most microorganisms were positively correlated to the contents of soil organic matter so that we might stabilize soil microbial communities and increase the abundance of beneficial bacteria by improving soil fertility. Notably, the research found that the abundance of Bacteroidetes increased in extremely harsh soils, and Bacteroidetes transferred and enriched a large amount of Ni ( Figure S2), which is consistent with Jiao et al. (2022) findings. Besides, Tang et al. (2019) found that Bacteroidetes also affected the availability of Cu and Zn. Okkeri and Haltia (1999; found that the zinctransporting P-type ATPase C was not only related to the transport of Zn but also to Cd and Pb, which is consistent with our genetic prediction. Similarly, nickel transport protein is a novel metal regulatory protein associated with heavy metal Ni transport (Dosanjh and Michel, 2006;Wang et al., 2014). In addition, there are functional genes secreting siderophore (Carroll and Moore, 2018), indoleacetic acid (Estenson et al., 2018) and ACC deaminase (Glick, 2005) in rhizosphere soil microorganisms. It is indicated that many plant growth- Schematic plot of the interaction between rhizosphere bacteria and plants in multi-HMs contaminated soil. Tryptophan synthase alpha/beta chain 6.71E-04a 5.82E-04b 4.89E-04c Nonribosomal peptide synthetases 8.82E-04b 7.33E-04b 1.21E-03a Nickel transport protein 1.72E-05a 3.72E-06b 3.47E-07b Manganese/Zinc-transporting P-type ATPase C 2.06E-05a 7.48E-06b 1.51E-05a promoting bacteria (PGPB) resist HMs in rhizosphere soil. They significantly improved plant growth in heavy metalcontaminated soils and could enhance heavy metal phytoremediation by binding super-enriched plants (Ahemad, 2019). Xu et al. (2022) found that Actinobacteriota and Gemmatimonadota promoted plant growth and fixed Cd in rhizosphere soil. Li et al. (2022) found that Sphingomonas promoted plant growth and degraded organometallic compounds to remediate soil HMs contamination. Similarly, Ham et al. (2022) found that Pseudopyroactor promoted plant growth and improved antioxidant capacity. PGPB significantly improve plant growth in heavy metal-contaminated soils (Ahemad, 2019) and may represent potential bacterial strain resources in plant-PGPR combined remediation of HMcontaminated soils. The presence of HMs significantly destabilizes network structure (Jiang et al., 2020;Caili et al., 2022). A more complex network represents that microbial activities and interactions are more active and intensive, which might have beneficial functions in the phytoremediation of HMs in soils (Hou et al., 2019). Plants consistently interact with a core set of microbes contributing to plant performance (Luo et al., 2022). We found that the network of LZ was the most complex ( Figure 3C), and Boehmeria nives (L.) showed higher TFs of all ten metals. Thus, a good rhizosphere ecological network might be helpful for plants to absorb heavy metals.

Conclusion
According to our investigation of BCF and TF for metals, LA preferred accumulating Cd, LZ preferred accumulating Cr and Sb, and LM preferred accumulating Cr and Ni. The bioaccumulating capacity for multi-metals of Artemisia argyi Levl was more potent, while the translocating capacity was weaker than that of Miscanthus floridulus (Lab.) and Boehmeria nives (L.). The dominant bacterial groups, ecological networks and soil properties (e.g., available phosphorus and moisture content) were different in the rhizosphere soil of these three plants. A higher proportion of positive links in the ecological network might enhance the metal uptake of plants. PICRUSt analysis and correlation tests indicated that genes related to the synthesis of proteins (e.g., manganese/zinc-transporting Ptype ATPase C, nickel transport protein and ACC deaminase) could also promote phytoremediation. The rhizosphere bacterial community contained genera, such as Truepera, Anderseniella, Paracoccus, Erythrobacter and Novosphingobium, which may represent potential bacterial strain resources for plant-microbes combined remediation of HM-contaminated soils.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.